#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include "HYPRE_struct_ls.h"
#include "_hypre_struct_ls.h"

#define MIN(a, b) (a) < (b) ? (a) : (b)
#ifndef TEST_CNT
#define TEST_CNT 1
#endif

int main(int argc, char * argv[])
{
    int my_pid, num_procs;
    HYPRE_Int glb_nrows;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);

    HYPRE_Int cnt = 1;
    const char * case_name = argv[cnt++];
    HYPRE_Int gx = atoi(argv[cnt++]), gz = atoi(argv[cnt++]);
    HYPRE_Int px = atoi(argv[cnt++]), pz = atoi(argv[cnt++]);
    const char * precond_name = argv[cnt++];
    if (strcmp(case_name, "OCEAN") == 0) {
        glb_nrows = gx * gz;
    }
    else if (strstr(case_name, "DEMO")) {
        glb_nrows = gx * gz;
    }
    else if (strstr(case_name, "upwind")) {
        glb_nrows = gx * gz;
    }
    else {
        if (my_pid == 0) printf("Error: Unrecognized case name %s\n", case_name);
        MPI_Abort(MPI_COMM_WORLD, -1);
    }
    HYPRE_Int lx, lz;
    const HYPRE_Int ndim = 2;

    MPI_Comm cart_comm;
    int cart_ids[ndim], pidx, pidy, pidz;
    HYPRE_Int ilower[ndim], iupper[ndim];// 注意i[2]在最外维，i[0]在最内维 所以令0: z, 1: x, 2: y
    {// Setup Cartesian comm
        HYPRE_Int dims[ndim] = {px, pz};// 注意这里保持一致！！！
        HYPRE_Int periods[ndim] = {0, 0};
        MPI_Cart_create(MPI_COMM_WORLD, ndim, dims, periods, 0, &cart_comm);
        MPI_Comm_rank(cart_comm, &my_pid);
        MPI_Comm_size(cart_comm, &num_procs);
        MPI_Cart_coords(cart_comm, my_pid, ndim, cart_ids);
        pidx = cart_ids[0];
        pidz = cart_ids[1];

        lx = gx / px;
        lz = gz / pz;
        ilower[0] = pidz * lz; iupper[0] = ilower[0] + lz - 1;
        ilower[1] = pidx * lx; iupper[1] = ilower[1] + lx - 1;
        int remain_x = gx - lx * px;
        int remain_z = gz - lz * pz;
        if (pidx < remain_x) {
            ilower[1] += pidx;
            iupper[1] += pidx + 1;
        } else {
            ilower[1] += remain_x;
            iupper[1] += remain_x;
        }
        if (pidz < remain_z) {
            ilower[0] += pidz;
            iupper[0] += pidz + 1;
        } else {
            ilower[0] += remain_z;
            iupper[0] += remain_z;
        }
        // 更新本进程的范围
        lx = iupper[1] + 1 - ilower[1];
        lz = iupper[0] + 1 - ilower[0];
        // printf("proc %d (%3d,%3d) (%3d,%3d)\n", my_pid, ilower[1], ilower[0], iupper[1], iupper[0]);
    }
    HYPRE_Init();
    HYPRE_StructGrid     grid;
    HYPRE_StructStencil  stencil;
    HYPRE_StructMatrix   A;
    HYPRE_StructVector   b, x, y;
    HYPRE_StructSolver   solver;
    HYPRE_StructSolver   precond;

    for (HYPRE_Int test = 0; test < TEST_CNT; test++) {
    HYPRE_Int num_diag = -1;
    {// Setup a grid
        HYPRE_StructGridCreate(cart_comm, ndim, &grid);// Create an empty 3D grid object
        HYPRE_StructGridSetExtents(grid, ilower, iupper);// Add a new box to the grid
        HYPRE_StructGridAssemble(grid);// a collective call finalizing the grid assembly.
    }
    if (strcmp(case_name, "OCEAN") == 0) {
        assert(sizeof(HYPRE_Real) == sizeof(double));
        num_diag = 5;
        HYPRE_Int offsets[num_diag][ndim] = {{0, -1}, {-1, 0}, {0,0}, {1, 0}, {0, 1}};
        {// Define the discretization stencil
            HYPRE_StructStencilCreate(ndim, num_diag, &stencil);
            // Define the geometry of the stencil
            for (HYPRE_Int entry = 0; entry < num_diag; entry++)
                HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
        }
    }
    else if (strstr(case_name, "DEMO")) {
        assert(sizeof(HYPRE_Real) == sizeof(double));
        if (strstr(case_name, "05")) {
            num_diag = 5;
            HYPRE_Int offsets[num_diag][ndim] = {{0, -1}, {-1, 0}, {0,0}, {1, 0}, {0, 1}};
            HYPRE_StructStencilCreate(ndim, num_diag, &stencil);// Define the discretization stencil
            // Define the geometry of the stencil
            for (HYPRE_Int entry = 0; entry < num_diag; entry++)
                HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
        }
        else if (strstr(case_name, "09")) {
            num_diag = 9;
            HYPRE_Int offsets[num_diag][ndim] = {{-1,-1}, {0, -1}, {1,-1}, {-1, 0}, {0,0}, {1, 0}, {-1,1}, {0, 1}, {1,1}};
            HYPRE_StructStencilCreate(ndim, num_diag, &stencil);// Define the discretization stencil
            // Define the geometry of the stencil
            for (HYPRE_Int entry = 0; entry < num_diag; entry++)
                HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
        }
    }
    else if (strcmp(case_name, "upwind") == 0) {
        assert(sizeof(HYPRE_Real) == sizeof(double));
        num_diag = 9;
        HYPRE_Int offsets[num_diag][ndim] = {{-1,-1}, {0, -1}, {1,-1}, {-1, 0}, {0,0}, {1, 0}, {-1,1}, {0, 1}, {1,1}};
        HYPRE_StructStencilCreate(ndim, num_diag, &stencil);// Define the discretization stencil
        // Define the geometry of the stencil
        for (HYPRE_Int entry = 0; entry < num_diag; entry++)
            HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
    }
    else if (strcmp(case_name, "upwind-2d5") == 0) {
        assert(sizeof(HYPRE_Real) == sizeof(double));
        num_diag = 5;
        HYPRE_Int offsets[num_diag][ndim] = {{0, -1}, {-1, 0}, {0,0}, {1, 0}, {0, 1}};
        HYPRE_StructStencilCreate(ndim, num_diag, &stencil);// Define the discretization stencil
        // Define the geometry of the stencil
        for (HYPRE_Int entry = 0; entry < num_diag; entry++)
            HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
    }

    assert(num_diag != -1);
    HYPRE_Int stencil_indices[num_diag];
    for (HYPRE_Int j = 0; j < num_diag; j++)
        stencil_indices[j] = j;
    {// Set up a Struct Matrix and Struct Vector
        HYPRE_Int tot_len = num_diag * lx * lz;
        HYPRE_Real * buf = (HYPRE_Real *) malloc (sizeof(HYPRE_Real) * tot_len);
        // 首先MPI-IO读入矩阵数据，或者自己构造数据
        if (strcmp(case_name, "OCEAN") == 0 || strcmp(case_name, "upwind") == 0 || strcmp(case_name, "upwind-2d5") == 0) {
            assert(sizeof(double) == sizeof(HYPRE_Real));
            MPI_File fh = MPI_FILE_NULL;// 文件句柄
            MPI_Datatype etype = (sizeof(HYPRE_Real) == 4) ? MPI_FLOAT : MPI_DOUBLE;
            MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
            HYPRE_Int size[3], subsize[3], start[3];
            size   [0] = gx       ; size   [1] = gz       ; size   [2] = num_diag;
            subsize[0] = lx       ; subsize[1] = lz       ; subsize[2] = num_diag;
            start  [0] = ilower[1]; start  [1] = ilower[0]; start  [2] = 0       ;// 注意这里要反过来

            MPI_Type_create_subarray(3, size, subsize, start, MPI_ORDER_C, etype, &read_type);
            MPI_Type_commit(&read_type);
            MPI_Offset displacement = 0;
            displacement *= sizeof(HYPRE_Real);// 位移要以字节为单位！
            MPI_Status status;

            int ret;
            ret = MPI_File_open(cart_comm, "mat.AOS.bin", MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
            if (ret != MPI_SUCCESS) {
                printf("Could not open file: ret %d\n", ret);
            }
            MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
            MPI_File_read_all(fh, buf, tot_len, etype, &status);
            MPI_File_close(&fh);
            MPI_Type_free(&read_type);

            double loc_dot = 0.0, glb_dot = 0.0;
            for (int i = 0; i < lx * lz; i++) {
                if (buf[i * num_diag + 2] == -1.0) {
                    continue;
                }
                // for (int d = 0; d < num_diag; d++)
                    loc_dot += buf[i * num_diag + 2] * buf[i * num_diag + 2];
            }
            MPI_Allreduce(& loc_dot, & glb_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
            if (my_pid == 0) printf("(A,A) = %.16e\n", glb_dot);
        }
        else if (strstr(case_name, "DEMO")) {
            const HYPRE_Int x_lbdr = (pidx == 0), x_ubdr = (pidx == px - 1),
                            z_lbdr = (pidz == 0), z_ubdr = (pidz == pz - 1);
            if (strstr(case_name, "05")) {
                for (HYPRE_Int i = 0; i < lx; i++)
                for (HYPRE_Int k = 0; k < lz; k++) {
                    HYPRE_Real * ptr = buf + (i * lz + k) * num_diag;
                    ptr[0] = (i == 0 && x_lbdr) ? 0 : -1.0;
                    ptr[1] = (k == 0 && z_lbdr) ? 0 : -1.0;
                    ptr[2] = 4.0;
                    ptr[3] = (k == lz-1 && z_ubdr) ? 0 : -1.0;
                    ptr[4] = (i == lx-1 && x_ubdr) ? 0 : -1.0;
                }
            }
            else if (strstr(case_name, "09")) {
                for (HYPRE_Int i = 0; i < lx; i++)
                for (HYPRE_Int k = 0; k < lz; k++) {
                    HYPRE_Real * ptr = buf + (i * lz + k) * num_diag;
                    ptr[0] = ((i == 0    && x_lbdr) || (k == 0    && z_lbdr)) ? 0 : -1.0;
                    ptr[1] =  (i == 0    && x_lbdr)                           ? 0 : -1.0;
                    ptr[2] = ((i == 0    && x_lbdr) || (k == lz-1 && z_ubdr)) ? 0 : -1.0;
                    ptr[3] =                           (k == 0    && z_lbdr)  ? 0 : -1.0;
                    ptr[4] = 8.0;               
                    ptr[5] =                           (k == lz-1 && z_ubdr)  ? 0 : -1.0;
                    ptr[6] = ((i == lx-1 && x_ubdr) || (k == 0    && z_lbdr)) ? 0 : -1.0;
                    ptr[7] =  (i == lx-1 && x_ubdr)                           ? 0 : -1.0;
                    ptr[8] = ((i == lx-1 && x_ubdr) || (k == lz-1 && z_ubdr)) ? 0 : -1.0;
                }
            }
            else MPI_Abort(MPI_COMM_WORLD, -509);
        }

        // 然后 传给结构化矩阵A
        HYPRE_StructMatrixCreate(cart_comm, grid, stencil, &A);// Create an empty matrix object
        HYPRE_StructMatrixInitialize(A);// Indicate that the matrix coefficients are ready to be set 
        HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, num_diag, stencil_indices, buf);
        HYPRE_StructMatrixAssemble(A);// a collective call finalizing the matrix assembly.
        
        free(buf);

        // 向量
        HYPRE_StructVectorCreate(cart_comm, grid, &b);// Create an empty vector object
        HYPRE_StructVectorCreate(cart_comm, grid, &x);
        HYPRE_StructVectorCreate(cart_comm, grid, &y);
        HYPRE_StructVectorInitialize(b);// Indicate that the vector coefficients are ready to be set
        HYPRE_StructVectorInitialize(x);
        HYPRE_StructVectorInitialize(y);

        tot_len = lx * lz;
        buf = (HYPRE_Real *) malloc (sizeof(HYPRE_Real) * tot_len);
        if (strcmp(case_name, "OCEAN") == 0 || strcmp(case_name, "upwind") == 0 || strcmp(case_name, "upwind-2d5") == 0) {
            assert(sizeof(double) == sizeof(HYPRE_Real));
            MPI_File fh = MPI_FILE_NULL;// 文件句柄
            MPI_Datatype etype = (sizeof(HYPRE_Real) == 4) ? MPI_FLOAT : MPI_DOUBLE;
            MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
            HYPRE_Int size[2], subsize[2], start[2];
            size   [0] = gx       ; size   [1] = gz       ;
            subsize[0] = lx       ; subsize[1] = lz       ;
            start  [0] = ilower[1]; start  [1] = ilower[0];// 注意这里要反过来

            MPI_Type_create_subarray(2, size, subsize, start, MPI_ORDER_C, etype, &read_type);
            MPI_Type_commit(&read_type);
            MPI_Offset displacement = 0;
            displacement *= sizeof(HYPRE_Real);// 位移要以字节为单位！
            MPI_Status status;

            int ret;
            ret = MPI_File_open(cart_comm, "b.AOS.bin", MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
            if (ret != MPI_SUCCESS) {
                printf("Could not open b file: ret %d\n", ret);
            }
            MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
            MPI_File_read_all(fh, buf, tot_len, etype, &status);
            MPI_File_close(&fh);
            HYPRE_StructVectorSetBoxValues(b, ilower, iupper, buf);

            ret = MPI_File_open(cart_comm, "x.AOS.bin", MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
            if (ret != MPI_SUCCESS) {
                printf("Could not open x file: ret %d\n", ret);
            }
            MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
            MPI_File_read_all(fh, buf, tot_len, etype, &status);
            MPI_File_close(&fh);
            MPI_Type_free(&read_type);
            HYPRE_StructVectorSetBoxValues(x, ilower, iupper, buf);
        }
        else if (strstr(case_name, "DEMO")) {
            for (HYPRE_Int i = 0; i < lx; i++)
            for (HYPRE_Int k = 0; k < lz; k++) buf[i * lz + k] = 1.0;
            HYPRE_StructVectorSetBoxValues(b, ilower, iupper, buf);
        
            for (HYPRE_Int i = 0; i < lx; i++)
            for (HYPRE_Int k = 0; k < lz; k++) buf[i * lz + k] = 0.0;
            HYPRE_StructVectorSetBoxValues(x, ilower, iupper, buf);
        }
        
        HYPRE_StructVectorAssemble(b);
        HYPRE_StructVectorAssemble(x);
        HYPRE_StructVectorAssemble(y);// 对y也要创建
        
        free(buf);
    }

    HYPRE_Real b_dot;
    b_dot = hypre_StructKrylovInnerProd(b, b);
    if (my_pid == 0) printf("(  b,   b) = %.20e\n",  b_dot);
    {// Check for DEBUG purpose
        HYPRE_Real x_dot, Ab_dot, Ax_dot;
        x_dot = hypre_StructKrylovInnerProd(x, x);
        if (my_pid == 0) printf("(  x,   x) = %.20e\n",  x_dot);
        HYPRE_StructMatrixMatvec(1.0, A, x, 0.0, y);
        Ax_dot= hypre_StructKrylovInnerProd(y, y);
        if (my_pid == 0) printf("(A*x, A*x) = %.20e\n", Ax_dot);
        HYPRE_StructMatrixMatvec(1.0, A, b, 0.0, y);
        Ab_dot= hypre_StructKrylovInnerProd(y, y);
        if (my_pid == 0) printf("(A*b, A*b) = %.20e\n", Ab_dot);
    }
    
    if (strcmp(case_name, "OCEAN") == 0 || strcmp(case_name, "upwind") == 0 || strcmp(case_name, "upwind-2d5") == 0) {
        HYPRE_StructGMRESCreate(cart_comm, &solver);
        HYPRE_StructGMRESSetKDim(solver, 10);
        HYPRE_StructGMRESSetMaxIter(solver, 1000);
        HYPRE_StructGMRESSetAbsoluteTol(solver, 0.0);
        // HYPRE_StructGMRESSetAbsoluteTol(solver, 1.0e-09);
        HYPRE_StructGMRESSetTol(solver, 1.0e-9);
        // HYPRE_StructGMRESSetTol(solver, 0.0);
        HYPRE_StructGMRESSetPrintLevel(solver, 2);
        HYPRE_StructGMRESSetLogging(solver, 1);

        if (strcmp(precond_name, "SMG") == 0) {
            if (my_pid == 0) printf("using SMG...\n");
            
            HYPRE_StructSMGCreate(cart_comm, &precond);
            HYPRE_StructSMGSetPrintLevel(precond, 10);
            HYPRE_StructSMGSetLogging(precond, 0);
            // HYPRE_StructSMGSetMemoryUse(precond, 0);
            HYPRE_StructSMGSetMaxIter(precond, 1);
            HYPRE_StructSMGSetTol(precond, 0.0);

            HYPRE_StructSMGSetZeroGuess(precond);
            // SMG的松弛类型是不能变的，它做的是一维半粗化，然后每个面做plane-solve，其中的plane-solve又来一遍半粗化和线松弛
            HYPRE_StructSMGSetNumPreRelax(precond, 1);
            HYPRE_StructSMGSetNumPostRelax(precond, 1);
            HYPRE_StructGMRESSetPrecond(solver,   HYPRE_StructSMGSolve,
                                                HYPRE_StructSMGSetup, precond);
        }
        else if (strcmp(precond_name, "PFMG") == 0) {
            if (my_pid == 0) printf("using PFMG...\n");
            HYPRE_StructPFMGCreate(cart_comm, &precond);
            HYPRE_StructPFMGSetPrintLevel(precond, 3);
            HYPRE_StructPFMGSetLogging(precond, 0);

            HYPRE_StructPFMGSetMaxIter(precond, 1);
            HYPRE_StructPFMGSetTol(precond, 0.0);
            HYPRE_StructPFMGSetZeroGuess(precond);
            HYPRE_StructPFMGSetRAPType(precond, 0);
            if (strstr(case_name, "upwind")) HYPRE_StructPFMGSetRelaxType(precond, 2);
            else                             HYPRE_StructPFMGSetRelaxType(precond, 1);
            HYPRE_StructPFMGSetNumPreRelax(precond, 1);
            HYPRE_StructPFMGSetNumPostRelax(precond, 1);
            HYPRE_StructGMRESSetPrecond(solver,   HYPRE_StructPFMGSolve,
                                                HYPRE_StructPFMGSetup, precond);
        }
        else if (strcmp(precond_name, "NOPC") == 0) { 
            if (my_pid == 0) printf("  NO precond!\n");
        }
        else {
            if (my_pid == 0) printf("  Error: Uncognized PC name: %s!\n", precond_name);
            MPI_Abort(MPI_COMM_WORLD, -101);
        }

        HYPRE_Real t_setup = MPI_Wtime();
        HYPRE_StructGMRESSetup(solver, A, b, x ); t_setup = MPI_Wtime() - t_setup;
        HYPRE_Real t_calc  = MPI_Wtime();
        HYPRE_StructGMRESSolve(solver, A, b, x);  t_calc  = MPI_Wtime() - t_calc;

        /* Get info and release memory */
        HYPRE_Int num_iterations;
        HYPRE_Real final_res_norm;
        HYPRE_StructGMRESGetNumIterations( solver, &num_iterations );
        HYPRE_StructGMRESGetFinalRelativeResidualNorm( solver, &final_res_norm );

        if (my_pid == 0) {
            printf("\n");
            printf("Iterations = %d\n", num_iterations);
            printf("Time cost %.5f %.5f %.5f %d\n", t_setup, t_calc, t_setup + t_calc, num_iterations);
            printf("Final Relative Residual Norm = %e\n", final_res_norm);
            printf("\n");
        }

        // Destroy solver and preconditioner
        HYPRE_StructGMRESDestroy(solver);
        if      (strcmp(precond_name, "SMG" ) == 0) HYPRE_StructSMGDestroy (precond);
        else if (strcmp(precond_name, "PFMG") == 0) HYPRE_StructPFMGDestroy(precond);
    }
    else if (strstr(case_name, "DEMO")) {
        HYPRE_StructPCGCreate(cart_comm, &solver);
        HYPRE_StructPCGSetMaxIter(solver, 1000);
        HYPRE_StructPCGSetTol(solver, 1.0e-09);

        HYPRE_StructPCGSetTwoNorm(solver, 1);
        HYPRE_StructPCGSetRelChange(solver, 0);
        HYPRE_StructPCGSetPrintLevel(solver, 2);
        HYPRE_StructPCGSetLogging(solver, 1);

        if (strcmp(precond_name, "SMG") == 0) {
            if (my_pid == 0) printf("using SMG...\n");
            HYPRE_StructSMGCreate(cart_comm, &precond);
            HYPRE_StructSMGSetPrintLevel(precond, 10);
            HYPRE_StructSMGSetLogging(precond, 10);
            HYPRE_StructSMGSetMaxIter(precond, 1);
            HYPRE_StructSMGSetTol(precond, 0.0);
            HYPRE_StructSMGSetZeroGuess(precond);
            // SMG的松弛类型是不能变的，它做的是一维半粗化，然后每个面做plane-solve，其中的plane-solve又来一遍半粗化和线松弛
            HYPRE_StructSMGSetNumPreRelax(precond, 1);
            HYPRE_StructSMGSetNumPostRelax(precond, 1);
            HYPRE_StructPCGSetPrecond(solver,   HYPRE_StructSMGSolve,
                                                HYPRE_StructSMGSetup, precond);
        }
        else if (strcmp(precond_name, "PFMG") == 0) {
            if (my_pid == 0) printf("using PFMG...\n");
            HYPRE_StructPFMGCreate(cart_comm, &precond);
            HYPRE_StructPFMGSetPrintLevel(precond, 3);
            HYPRE_StructPFMGSetLogging(precond, 0);
            HYPRE_StructPFMGSetMaxIter(precond, 1);
            HYPRE_StructPFMGSetTol(precond, 0.0);
            HYPRE_StructPFMGSetZeroGuess(precond);
            // HYPRE_StructPFMGSetRAPType(precond, 1);
            HYPRE_StructPFMGSetRelaxType(precond, 1);// 0: 655次，1: 70次（如果搭配RAPType=1则要102次），2：93次，3：89次（2和3不能使用RAPType=0）
            HYPRE_StructPFMGSetNumPreRelax(precond, 1);
            HYPRE_StructPFMGSetNumPostRelax(precond, 1);
            HYPRE_StructPCGSetPrecond(solver,   HYPRE_StructPFMGSolve,
                                                HYPRE_StructPFMGSetup, precond);
        }
        else if (strcmp(precond_name, "NOPC") == 0) { 
            if (my_pid == 0) printf("  NO precond!\n");
        }
        else {
            if (my_pid == 0) printf("  Error: Uncognized PC name: %s!\n", precond_name);
            MPI_Abort(MPI_COMM_WORLD, -101);
        }

        HYPRE_Real t_setup = MPI_Wtime();
        HYPRE_StructPCGSetup(solver, A, b, x ); t_setup = MPI_Wtime() - t_setup;
        HYPRE_Real t_calc  = MPI_Wtime();
        HYPRE_StructPCGSolve(solver, A, b, x);  t_calc  = MPI_Wtime() - t_calc;
        /* Get info and release memory */
        HYPRE_Int num_iterations;
        HYPRE_Real final_res_norm;
        HYPRE_StructPCGGetNumIterations( solver, &num_iterations );
        HYPRE_StructPCGGetFinalRelativeResidualNorm( solver, &final_res_norm );
        if (my_pid == 0) {
            printf("\n");
            printf("Iterations = %d\n", num_iterations);
            printf("Time cost %.5f %.5f %.5f %d\n", t_setup, t_calc, t_setup + t_calc, num_iterations);
            printf("Final Relative Residual Norm = %e\n", final_res_norm);
            printf("\n");
        }
        // Destroy solver and preconditioner
        HYPRE_StructPCGDestroy(solver);
        if      (strcmp(precond_name, "SMG" ) == 0) HYPRE_StructSMGDestroy (precond);
        else if (strcmp(precond_name, "PFMG") == 0) HYPRE_StructPFMGDestroy(precond);
    }

    // 计算真实残差
    HYPRE_StructVectorCopy(b, y);// y = b
    HYPRE_StructMatrixMatvec(-1.0, A, x, 1.0, y);// y += -A*x
    HYPRE_Real r_dot;
    r_dot = hypre_StructKrylovInnerProd(y, y);
    if (my_pid == 0) printf("\033[1;35mtrue ||r||/||b||= %20.16e\033[0m\n", sqrt(r_dot/b_dot));

    /* Free memory */
    HYPRE_StructGridDestroy(grid);
    HYPRE_StructStencilDestroy(stencil);
    HYPRE_StructMatrixDestroy(A);
    HYPRE_StructVectorDestroy(b); HYPRE_StructVectorDestroy(x); HYPRE_StructVectorDestroy(y);
    }// test loop

    /* Finalize HYPRE */
    HYPRE_Finalize();

    /* Finalize MPI */
    MPI_Finalize();
    return 0;
}